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Abstract. Thanks to PfafSan techniques, we study the Renyi entanglement entropies 

and the entanglement spectrum of large subsystems for two-dimensional Rokhsar- 
Kivelson wave functions constructed from a dimer model on the triangular lattice. By 
including a fugacity t on some suitable bonds, one interpolates between the triangular 
lattice {t = 1) and the square lattice {t = 0). The wave function is known to be a 
massive Z2 topological liquid for t > whereas it is a gapless critical state at i = 0. 
We mainly consider two geometries for the subsystem: that of a semi-infinite cylinder, 
and the disk-like setup proposed by Kitaev and Preskill [Phys. Rev. Lett. 96, 110404 
(2006)]. In the cylinder case, the entropies contain an extensive term - proportional to 
the length of the boundary - and a universal sub-leading constant s„ (t) . Fitting these 
cylinder data (up to a perimeter of i = 32 sites) provides s„ with a very high numerical 
accuracy (10~^ at f = 1 and 10~^ a,t t = 0.5). In the topological Z2 liquid phase we 
find Sn{t > 0) = — ln2, independent of the fugacity t and the Renyi parameter n. At 
<: = we recover a previously known result, s„(t = 0) = — | ln(n)/(n — 1) for n < 1 
and Sn{t = 0) = — ln(2)/(n — 1) for n > 1. In the disk-like geometry - designed to 
get rid of the boundary contributions - we find an entropy (i > 0) = — In 2 in the 
whole massive phase whatever n > 0, in agreement with the result of Flammia et al. 
[Phys. Rev. Lett. 103, 261601 (2009)]. Some results for the gapless limit R^it ^ 0) 
are discussed. 
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1. Introduction 

It is now widely recognized that the entanglement entropy is a useful quantity to probe 
many-body quantum states. It can be used to detect critical states in one-dimensional 
chains, through the celebrated logarithmic divergence [H [2l [3] Hj. In two dimensions it 
can be a used to characterize (massive) topologically ordered states. In particular, it 
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allows to distinguish a topological wave function from a more conventional disordered 
and featureless state. In a gapped phase the entanglement entropy of a large subsystem 
contains a contribution proportional to the length (in two dimensions) of its boundary 
plus a subleading term S'topo which contains some information about the nature of the 
phase. In a state with topological order, this subleading term is related to the total 
quantum dimension, that is to the content in elementary excitation O |6l [7] . This idea 
has been successfully applied to some fractional quantum hall states [HllHlIin]- Extracting 
the subleading term in lattice models is not a trivial task [HI [7] but it was first shown 
to be feasible using quantum dimer wave functions on the triangular lattice [illj. Since 
the work of Moessner and Sondhi jjl] this type of states have been intensively studied 
since they offer some rather simple realization of topologically ordered states with non 
trivial finite-size effects and finite correlation length (contrary to toric-code like models 

mm)- 

In this paper we also consider some dimer wave functions - named after Rokhsar and 
Kivelson (RK) [15] - which are linear superposition of fully packed dimer coverings on 
the triangular lattice. By including a fugacity on some suitable bonds, one continuously 
interpolates between the triangular lattice [t = 1) and the square lattice {t = 0). In 
the triangular case the wave function is known to be a massive Z2 topological liquid 
fn\ Uni [IZ] whereas it is a gapless critical state at t = [IS]. Exploiting previous results 
[TTl [18] on the reduced density matrix (RDM) of RK states, we can obtain not only 
the entanglement entropy but also the full entanglement spectrum on large systems. 
Using extensively the Pfaffian formulation of the classical dimer partition function [19] , 
as well as some perturbation theory for determinant [201 HB] we perform calculations in 
the thermodynamic limit while keeping the boundary length finite. 

In the cylinder geometry we can treat the infinite height limit and perimeters up to 
L = 32 (38 at t = 0). In the disk-like geometry proposed by Kitaev and Preskill [6J, we 
perform exact calculation for disks of radii up to p ~ 4.5 lattice spacings embedded in a 
infinite system, therefore extending significantly the previous entanglement calculations 
on triangular dimer wave functions [TT]. This technique allows to confirm the value 
•S'topo = — lii(2) with high precision in the whole massive phase (not only at the triangular 
point t = 1). This value turns out not to depend on the Renyi parameter, in agreement 
with the argument by Flammia et al. [21j. We also discuss the structure of the 
entanglement spectrum, showing that it contains a non-degenerate "ground-state" and 



a gap. In Sec. |4.6[ a micro-canonical point of view is used to relate the density of states 
of the entanglement spectrum and the Renyi entanglement entropies. 

When t = the dimers are restricted to the bonds of a square lattice. Although 
non-genericjlj such critical RK wave-functions associated to some conformally invariant 
critical points are useful since they offer one of the very few situations where one 
can study the entanglement in a critical wave-functions in more than one dimension 
[25] |26| [T8l [28] [29] . Another point of view is that, for long cylinder geometries, the 



I They correspond to fine tuned multi-critical points 
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entanglement in these two-dimensional systems is related to the Shannon entropies in - 
now generic - quantum critical chains [HI [30l [311 132] . The sub- leading constant in the 
cylinder geometry depends on the compactification radius |26l [THl EH |29] and shows a 
singularity at some critical value of the Renyi parameter [32]. The result in a Kitaev- 
Preskill geometry is less clear and we discuss our numerical results at the end of Sec. 

El 

2. Entanglement entropy as a Shannon entropy 

After a brief introduction to dimer RK wave-functions [15], we review how one can 
construct the RDM and Schmidt decomposition for these states. 

2.1. Rokhsar-Kivelson wave functions 

We start from a classical two-dimensional hard-core dimer model on a triangular lattice, 
with fugacity t on "diagonal" links (Fig. [T]). This fugacity allows to interpolate between 
the square lattice (t = 0) and the isotropic triangular lattice {t = 1). 



Figure 1. Triangular lattice with cylindrical boundary conditions {Lx — Q, Ly = 5). 
Each "diagonal link" (dotted lines) has fugacity t, the others have fugacity 1. 

The classical partition function of this system reads 

£ = e"'^*-'^^ = diagonal dimers ^-|^^ 

c c 

where the sum runs over all dimer coverings c. When t = (square lattice), the model 
is known to be critical [201 133] , its long range behavior is described by a compact free 
field [211 [35]. Otherwise it has a finite correlation length [121 [161 [17] . An Hilbert space is 
then constructed by associating a basis state |c) to each classical dimer configuration c. 
Different classical configurations correspond to orthogonal states. The RK wave function 
is the normalized linear combination of all basis states with an amplitude equal to the 
square root of the classical weight : 

|/2K) = -^$^e-^(^)/^|c). (2) 

* c 

Following Henley [37] one can construct some local Hamiltonians for which Eq. |2] 
is an exact ground state, but the precise form of these Hamiltonians will not be used in 
the following. 
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2.2. Renyi entanglement entropy 

We divide the system into two parts A and B. Each subsystem is a set of bonds, and its 
degrees of freedom are the corresponding dimer occupancies. The RDM of A is obtained 
by tracing over the degrees of freedom in B: 

Pa = Ttb\RK){RK\ (3) 

Then, the Renyi entanglement entropy is defined as 

5„ = ^lnTrp^, (4) 

where n is not necessarily an integer. Two limits are of interest. For n 1, Sn reduces 
to the Von Neumann entanglement entropy : 

5i = 5^^ = -Trp^lnpA (5) 

For n — )■ oo, only the largest eigenvalue Pmax of the RDM matters : 

Soo = - Inpmax- (6) 

This quantity is also called single copy entanglement. To compute all the Renyi 
entropies, we need all the eigenvalues of the RDM. In the following, we shall see that 
calculating each eigenvalue amounts to solving a combinatorial problem. The procedure 
has been discussed in details elsewhere [HI HH] and is recalled below for completeness. 

2.3. Schmidt decomposition 

We consider the geometry of an infinite cylinder cut into two parts, as in the left of 
Fig. [2j The reasoning is the same for the other geometries we considered. The sites 
which touch a bond in A and an bond in B (red circles in Fig. |2]) are called boundary 
sites. We assign a spin aj to each boundary site : aj =t if the site is occupied by a 
dimer in A, aj =4, if it is occupied by a dimer in B. We denote by 

\i) = \cTi,a2,. . . ,aLj (7) 

the whole spin configuration at the boundary. 

Now, let Sf^ (resp. S^) be the set of dimer configurations in A (resp. B) compatible 
with \i) at the boundary. Thanks to the hardcore constraint, they share no common 
element : 

Sfn£,'^ = (l} (8) 
Each configuration c can be written as 

c = aUb ,aeSf ,beSf (9) 
and the energy decomposed as 

E{c) = EA{a) + Esib) (10) 
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Figure 2. (Color online) Partition of the lattice in two subsystems A (red bonds) and 
B (blue bonds). Left: the subsystems A and B are semi-infinite cylinders. Boundary 
sites are marked by filled red circles. Each boundary site can either be occupied by a 
dimer in A (spin f), or a dimer in B (spin |). 



This allows to write the RK state as : 



E 



'EA{a)/2\ 



a) 



X 



E 



-Ss(b)/2 



Defining a new normalized set of RK states in A and B 
|RKf) = -^5:e-i--(")|a), 



with 



E 



Eq. (|2]) becomes 

\RK) = Y,^^RKt)\RKf), 



with 



Vi 



z ■ 



(11) 

(12) 

(13) 
(14) 

(15) 
(16) 



Eq. 15 is actually the Schmidt decomposition of the RK state (the orthogonality of the 



Schmidt vectors is guarantied by Eq. |8]), and the {p,} are the eigenvalues of the RDM: 

PA = Y.V^\RKt)mtV (17) 

i 

This way, one can obtain the Renyi entropy : 



s~~>"(Ep.") 



The entanglement entropy calculation has been reduced to finding some probabilities 
in the classical dimer problem. In the next section we will show that, using standard 
Pfaffian techniques, one can obtain exact formulae for the pj. 



CONTENTS 



7 



3. Classical probabilities 

3.1. Pfaffian 

The Pfaffian of a {2n x 2n) antisymmetric matrix M is defined as 

PiM= J2 e(7r)M,ri7r2M,r37r4 • • • Mt,,^_,t,,^, (19) 

where e{n) denotes the signature of a permutation vr. The sum runs over all permutations 
of {1,2,..., 2n} satisfying the constraints 

T!'2i-i < 7r2i , 1 <i <n ^2Q^ 
7r2i-i < vr2i+i , 1 < i < n - 1 
A very important relation is 

(PfM)^ = detM, (21) 

It is especially useful because it allows to compute the Pfaffian numerically in a time 
proportional to using standard determinant routines (and sometimes analytically). 

3.2. Kasteleyn theory 

The problem of enumerating dimer configurations on a planar lattice is a classic 
combinatorial problem, which was solved independently by Kasteleyn [12] and 
Temperley and Fisher |38]. We consider the case t = 1 for simplicity but the 
generalization to any t is straightforward. For any planar graph, the partition function 
(number of dimer coverings) is given by 

Z=|Pf/C|, (22) 

where /C is an antisymmetric matrix constructed as follows. Putting arrows on all the 
links, a matrix element of JC is 

{+1 if the arrow points from i to j 
— 1 if the arrow points from j to i (23) 
if i and j are not nearest neighbors 

The Kasteleyn matrix must also satisfy the clockwise- odd rule : the product of the 
arrows orientations (±1) around any elementary plaquette (running clockwise) has to 
be —1. Kasteleyn showed that i) such a matrix /C exists for any planar graph and ii) it 
insures that all terms in the sum have the same sign (the signature of the permutation 
always compensate that of the product of matrix elements). It is immediate to check 
that ii) implies Eq 



22 



A Kasteleyn matrix obeying Eq. |22] can also be found for cylindrical boundary 
conditions. An example for the triangular lattice with cylindrical boundary condition^ 
is shown in Fig. [3j 

§ In the case of toroidal boundary condition the situation is shghtly more compHcated, and the number 
of dimer covering is given by a Hnear combination of four Pfaffians, see Ref. 39 for more details. 
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In the following we will demonstrate how each probability pi can be computed as 
a determinant, taking the example of the cylinder geometry. 




Figure 3. Kasteleyn orientation of the [L^ = 6,Ly = 5) lattice (a weight t is given 
to "diagonal" links). Blue arrows: orientation of the bonds. Green : bonds present 
because of periodic boundary conditions along the x— axis (see Ref. |3£|). Their 
orientations are reversed compared to their "bulk" counterparts. 



3.3. Classical probabilities 



To find the probabilities of Eq. 16, we need to compute Zf^Z^, which is the 
partition function restricted to dimer configurations compatible with the boundary 
spin configuration \i) = \ai, . . . ai^) ■ It can be evaluated as the Pfaffian of a modified 
Kasteleyn matrix 

Zfzf = Pf (24) 

where /C*^*^ is deduced from /C by removing the appropriate links in a simple way. If 
aj =t, a dimer emanating from the boundary site j has to be in A, and we remove 
links in B emanating from site j. If aj =1 we remove links in A emanating from site j. 
See Fig. |4]for two examples, one with the boundary configuration \i) = \ tiitti) and 
one with \i) = Itttttt)- The computation of any such probability apparently requires 
the ratio of two L^Ly x L^Ly determinants. However, using a known trick [2D], the 
computation can be greatly simplified. 

3.4. Perturbation theory for determinants in an infinite system 

Following Ref. [20], pf may be written as 

p2 = det(l + JC-^S'^^) , = /C(^) - /C (25) 

The important point is that the matrix element S^'^J, is non zero only if the link r o r' 
has been removed. Then, a matrix element of }C~^£^^^ is 

(^~'^'^)rr' = E^r-X^:'- (26) 
s 

It is non-zero only if r' is a site belonging to a removed link. We name these sites "vicinity 
sites", and they of course depend on the boundary configuration \i). A boundary site is 
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Figure 4. Two examples with |i) = Itlltti) on the left, and \i) = Itttttt) on the 
right. Filled red circle: boundary site occupied by a dimer in A (spin ^). Empty red 
circle: boundary site occupied by a dimer in B (spin J,). To ensure that a boundary site 
be occupied by a dimer in A (resp. B), all edges in B (resp. A) coming from this site 
have to be removed. Notice that after the removal, A and B are disconnected, black 
circles filled in grey are sites which are connected to a boundary site through a link 
that has been removed. As explained in the text, the size of the determinant is given 
by the number of circles, pi is therefore a 16 x 16 determinant for the configuration 
on the left and a 12 x 12 determinant for the configuration on the right. 



automatically a vicinity site, but the converse is not true however. If we denote by Ei the 
set of vicinity sites and by rii their number, /C~^£*-*^ is a L^^Ly x Lr^Ly matrix, but only 
Ui columns are non identically zero. Then, using the antisymmetry of the determinant, 
any cell with indices r and r' not both in E^ can be set to zero by appropriate linear 
combinations of rows and columns. Therefore, the determinant may be computed as its 
restriction to the sites in E^. 

P? = det((l + /C-i^«)l^j (27) 

This so called "perturbation theory for determinants" has been previously used in 
Ref. [20] to compute exactly the monomer-monomer correlation on the square lattice 
in the thermodynamic limit {L,Ly — )■ oo), and further extended in Ref. ^6] to the 
triangular lattice. For computational purpose this is a huge simplification, because the 
size of the determinant has been reduced from L^Ly to rii ~ 0{Lx), and the total 
system we are interested in can be infinite {Ly — )■ oo). Contrary also to the transfer 
matrix approach [18] , this method allows us to treat any shape of boundary. This will 
be particularly useful while studying the geometry proposed by Kitaev and Preskill [6] . 

For this to work we also need to compute exactly certain matrix elements of the 
inverse Kasteleyn matrix /C^^. This can be done using standard Fourier and integral 
techniques, see [Appendix A 



Let us now specify the case of the (infinitely long) cylinder geometry cut into 
two parts. An example of spin configuration is shown in Fig. |4| where boundary sites 
are represented by red circles (filled or empty depending on the spin). Other vicinity 
sites are circles filled in grey. It is easy to check that 2Lx < Ui < SL^ for all boundary 
configurations. Since there are a priori 2^"= boundary configurations and each probability 
is of complexity ~ n?, the Renyi entropy can be evaluated in a time ~ x 2^^. This 
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Figure 5. Number of correct digits in the numerical estimate of the topological 
constant, as a function of the number of boundary sites. For the cylinder geometry 
we show the data for t = 1 (red circles) and t — 0.5 (blue triangles). The number of 
boundary sites is just L in this case, and the estimate is obtained by a fit to ai + si for 
two even consecutive values of L. The convergence to the correct value is exponentially 
fast, with an effective correlation length close to the dimer-dimer correlation length 
(which can have an imaginary part [161 117) , hence the oscillations we observe). For 
comparison we also show the data in the Kitaev-Preskill geometry, slightly anticipating 
on Sec. |6l 



allow US to go to relatively large system sizes of order ~ 30. 
4. Results for the infinite cylinder 

When the height Ly is infinite, the entropies S'„ only depend on the perimeter = L. 
As usual, the leading term is non universal and scales with L, and we are interested in 
the first subleading contribution s„: 



4.1. Topological entanglement entropy and Renyi index 

For gapped topological wave functions, the subleading constant si in the Von Neumann 
entropy has been shown to be related to the content of the phase in terms of 
fractionalized particles, and to the total quantum dimension D in particular [HI [?]: 
Si = —ln{D). In the original works the subleading constant si was extracted by 
combining the entropies of different subsystems in a planar geometry. We show here 
that the subleading term can be extracted in a - somewhat simpler - cylinder geometry 
(see also [TU]). 
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For t > the present dimer wave-functions realize the simplest topological phase, 
the so-called Z2 liquid with quantum dimension D = 2. One therefore expects to have 
Si = — In 2 in the whole topological phase. So far, this has only been checked numerically 
at t = 1 [iTj- In addition, Ref. [21] argues that this topological entanglement entropy 
is independent of the Renyi index n. We present here some results for infinitely high 
cylinders for various values of t and n, which support this result. The convergence 
to the topological entropy is exponentially fast, as can be seen in Fig. |5} For generic 
values of t and n, this allows us to get this constant with a very high accuracy: for 
example at t = 1 our best estimate is \si{t = 1) + ln2| ~ 10^^. It is widely believed 
that in massive phases the topological entropies (subleading terms) are independent of 
short-range correlations, but this is not proven. The present results, which strongly 
indicate that s„ = —ln{2) for any t > 0, therefore brings some additional support to 
the robustness of topological entropies. In general finite-size effects get larger when 
increasing n at fixed t, and it is more advisable to numerically study low-n Renyi 



entropies. However, as is shown in [Appendix B.2[ the calculation for n — > 00 simplifies 
greatly, and the result Soo(^ > 0) = — ln2 can even be obtained rigorously. We further 
discuss this result in Sec. |4.5[ At fixed n the convergence is also less clear when t 
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Figure 6. Sub- leading constants s„(t) for 3 different values of the Renyi parameter 
(n = 0.5, 1, 1.5). For each t and n, s„(t) is extracted from S„{L) using two consecutive 
even values of L (up to L = 32). In the thermodynamic limit the results are expected 
to converge to s„(f) = — ln2 for all n > and t > 0. 

is small since the correlation length ^(t) diverges when approaching t = and the 
finite-size effects become very important when L > ^(t). Still, the curve s„(t > 0) 
approaches — In 2 when L — )■ 00. The data plotted in Fig. [6] are indeed compatible with 
Sn(t) = — ln2 for all n = 0.5, 1, 2 and t > 0. The scaling close to t = will be discussed 
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later in Sec. 14.31 

4-2. Thermodynamical entropy 

The behavior for large values of the Renyi index n is displayed in Fig. [s] (triangular 
dots). Although it is roughly constant and close to — ln(2), due to the finite-size of 
the system there are some visible deviations for n > 3. This is even more visible if we 
consider a slightly different entropy, 5^, defined as: 

S^= (l-a„)ln(Z„) (29) 
Zn=J2P^ (30) 

i 

which can also be written as the Shannon entropy associated to the normalized 
probabilities pf. 

S^ = -J2P^HP.) With P.= ^. (31) 

Both entropies match at n = 1 (Sj^^^ = Sn=i) and are simply related otherwise: 
5^ = (1 — ndn){{l — n)Sn)- The "thermodynamic" entropy 5*"^ has also a leading 
term 0{L) and a sub-leading term, s^. The extensive (and non-universal) part is 
plotted in Fig. [7]as a function of the "temperature" 1/n. To stress the similarity with 
usual statistical mechanics, we also plotted the associated "specific heat" defined as a 
derivative of S'^: = —n^^. 

The sub-leading term is plotted in Fig. [s] (crosses). It is very close to — ln(2) 
at small n, but goes to = when n — > oo. This is indeed expected since the 
thermodynamic entropy S'ij^^^ - which corresponds to zero "temperature" - is equal 
to the log of the degeneracy of the configuration with the highest probability, which 
is non-degenerate in our case. However, the crossover from — ln(2) to takes place at 
values of n which are larger and larger when L — )■ cxd. This can be checked in the inset 
of Fig. [8| where the numerical data appear to be correctly fitted by 

^n»ina)~^'exp(-nA) (32) 



'n>ln(L) 
'n<ln(L) 

where A ~ 1.32 is the entanglement gap at t = 1. We finally note that the calculation 



Sn«in(L) ~ - ln(2) (33) 



of Pmax given in Sec. Appendix B.2 proves rigorously that limL_>oo linin^oo ■§„ = — ln(2) 



and lim^^oo linin^oo 3^ = 0. 

4-3. Scaling when t — > and L — > oo with fixed L ■ t 

The critical point t = has already been studied |18l [32] and is known to give: 



, \nR- , < n < 1 

n—1 ' 
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Figure 7. Thermodynamic entropy per site S^/L (monotonously increasing, right 
axis) and its associated "specific heat" (peaked at n ~ 0.25, left axis) C„ = —ndS^/dn. 
Fugacity t = 1. 



S 



;(o)H . (35) 



where the compactification radius is i? = 1 (free fermions) for the present dimer wave- 
functions, but could be tuned by adding some dimer-dimer interactions [55] . 

The correlation length ^(t) diverges as ^(t) ~ when t 1 p5]. In Fig. |9] 
we plot the subleading constant Sn(t,L) as a function oi L ■ t L/^(t). It appears 
that, for a given value of n, the data curves corresponding to different values of t and 
L approximately collapse onto each other. This shows that, when the system size L 
is much bigger than the correlation length C,{t) ~ t~^, we find the correct topological 
entanglement entropy s„ = — ln(2). On the other hand, when L is of the same order of 
magnitude than ^(t) (and much larger than the lattice spacing) s„ turns out to be some 
non-trivial function of n and L ■ t. When L ■ t — )■ the system effectively behaves as a 



critical system of dimers on a square lattice and s„ converges to Eq. 34, as expected 



4- 4- Entropy of a zig-zag line 



As explained in Sec. 2^, the eigenvalues of the RDM of a half infinite cylinder are the 
classical probabilities of the "spin" configurations \i) = \ai,a2, ■ ■ ■ ,o'i,). But one may 
also consider a zig-zag line and the probabilities pa of the dimer configurations on that 
lines. The "spins" are now replaced by the dimer occupancies (say or 1) of the zig-zag 
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Figure 8. Large n behavior of the subleading constant s„(t = 1) of the Renyi entropy, 
and s^{t = 1), the subleading constant of the thermodynamical entropy. They both 
give — ln(2) for smaU n, but differ for large n. This is a finite-size-effect: as shown 
in the inset, ~ exp(— An) for large n. We thus have s c:^ =~ — ln(2) up to 
n - ln(L). 



bonds. Theses probabilities can be computed using exactly the same perturbed-Pfaffian 
method as before. However, in terms of entanglement, the entropy we compute is that 



of a the "zig-zag" chain shown in the right of Fig. 10 Although the probabilities are 
computed in a very similar way, this calculation does not describe the entanglement of 
a two-dimensional subsystem, but that of a one-dimensional line winding around the 
cylinder. 

The associated entropies, already considered in Ref.[TT], have a leading term 
proportional to L and a subleading contribution of order The results, plotted 

in Fig. [T0| show that the subleading constant si has a dependence on t and system size 
L which is very similar to that of the half-cylinder entropy. It is possible that, as a 
function of L ■ t, the zig-zag line and half-infinite cylinder converge to the same curves 
for sufficiently large L. In any case, the zig-zag results clearly converges to — ln(2) in 
the thermodynamic limit for t > 0. 

One may ask if the zig-zag entropy would also give access to the quantum dimension 
for a general topologically ordered wave-function (not of RK type, and even not based 
on dimers). We believe that it is not the case. The present dimer RK states enjoy 
a special property: once the dimer occupancies are fixed along the zig-zag chain, the 
upper and lower half-cylinders are completely decoupled. For this reason, the entropy 




-1.2 I ' ' ' ' 1 

2 4 6 8 10 



L.t 



Figure 9. Sub-leading constants s„, as a function of t x L. For each values of n, 
the data corresponding to different values of t and L appear to be well described by a 
function oi t x L only. 

of the zig-zag chain is very close to that of a half cylinder. This would not hold for more 
generic states and a thick strip (sufficiently large compared to the correlation length) 
would be probably required to access the quantum dimension in general. 




2 4 6 8 10 



L.t 



Figure 10. Sub-leading constants for the entanglement entropy calculated numerically 
in two geometries : half-infinite cylinder and zig-zag strip (see text). 
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4-.5. Infinite Renyi and bipartite fidelity 

As already emphasized, the infinite-n Renyi hmit selects the largest eigenvalue of the 
RDM, which is the probability of the most likely configuration in the dimer language: 

5'oo = -lnpmax (36) 

For the cylinder geometry the corresponding boundary configuration |Vax) is 
particularly simple (see Fig. |4]for a graphical representation): 

|Vax) = |tt •••!), (37) 
and Pmax can be expressed as a ratio of simple partition functions: 

P.„ = Mm lY'-.f'/f , (38) 

where Zcyi{L, h) is the partition function for dimers on a finite cylinder of length L and 
height h. As detailed in Appendix B , we then find the following expression for 5*00 

l<m<L/2 , . 2 , , \ 

m 1 / I 1 sm ft — tcos/c \ 

^ \ 2 2y/t^ + sm^ k + sm^ k ' 



(2m — l)-K 



L 



from which one can extract the sub-leading constant 

Soo{t) = \ ? ' (40) 
^ ' 1 -ln2 , t > 0. ^ ' 

This result has already been mentioned in Sec. |4.1 The entropy 6*00 can also be 
considered from a different point of view. \RK) is the ground state of the Rokhsar- 
Kivelson Hamiltonian, and lives on a cylinder of length L and height h. This 
Hamiltonian may be written as 

H = H^uB = Ha + Hb + H^^'l (41) 

where Ha (resp. Hb) is the Rokhsar-Kivelson Hamiltonian restricted to sites in A (resp. 
B). We have [Ha, Hb] = and H^^^ contains all the interactions between A and B. If 
we denote by \A) (resp. \B)) the ground-state of Ha (resp. Hb), \A® B) = \A) (g) \B) 
the ground-state of Ha + Hb and by \AU B) = \RK) the ground state of Haub, then 
Pmax can be reformulated as 

Pmax =\{AU1 

Taking minus the logarithm we get 



Pm..= \{AUB\A0B)f (42) 



S^ = -\n\{AUB\A0B)f (43) 



The r.h.s of Eq. |43jhas been studied in Ref. [40] under the name logarithmic bipartite 
fidelity (LBF). The (infinite) Renyi entanglement entropy and the LBF are a priori not 
related, but we find that they are simply equal for this particular RK wave-function. In 
other words, performing a Schmidt decomposition on the total wave function \AU B), 
the Schmidt state with the highest Schmidt value is nothing but the ground state of 
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Ha + Hbi the RK Hamiltonian where all interactions between A and B were switched 
off. 

However, this relation does not hold exactly in general. For instance, in the Kitaev- 
Preskill or Levin- Wen geometry the boundaries are not straight and in that case the 
boundary dimer configuration | Vax) is not as simple as for the cylinder. Still, as pointed 
out in the equivalence between the LBF and 6*00 can hold for some more complex 
topological states such as the string nets states constructed by Levin and Wen We 
expect that for a generic {i.e. non RK) gapped state, the sub-leading term in the LBF 
and 5*00 should be the same in the thermodynamic limit (although, due to some mismatch 
at short distances, the extensive terms will differ). The argument is as follows: starting 
from a string net state where the correspondence works, we adiabatically modify the 
wave function toward the state we are interested in (without closing the gap). Doing 
so it is natural to expect that only the short-distance properties of the entanglement 
will be modified (hence the ~ L term) but not the sub-leading constant Sqo which is 
expected to be free from the contribution of local correlations. Although the robustness 
to changes in local correlations is is not proven in general, we provide in [Appendix B| a 



rigorous proof that the subleading term s^a is equal to — In 2 in the whole massive phase 
of the model (t > 0). 



4.6. Entanglement gap and entanglement spectrum 




Figure 11. Entanglement spectrum for L = 12 and for t = 1,0.7,0.3,0 from left to 
right. 

The spectrum of the RDM contains some rich information about the system. 
Looking at such spectra has been particularly fruitful in the context of the quantum 
Hall effect (QHE), where the entanglement spectrum was shown [H] to reflect some 
properties of the chiral gapless excitations which can propagate along an edge [12]. 
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With RK wave functions the RDM eigenvalues are simple classical probabilities and we 
thus have a relatively easy access to the entanglement spectra of large systems. 



2.4 




1 2 3 4 5 



t 



Figure 12. Entanglement gap as a function of t. It is maximum at t = (square 
lattice) a decreases slowly to zero when t — > oo. Except very close to t = (inset) the 
curves for L = 16 and 20 are practically indistinguishable at the scale of the figure, 
signaling negligible finite-size-effects. 

Such spectra shown in Figs. [TIHT2| where the probabilities Pi have been converted 
to "energies": Ei = — ln(pj/pmax)- The first observation is that these spectra have a 
unique ground-state and a gap A = to the first "excitation" . This is true not only 
in the Z2 liquid {t > 0) but also for the critical RK wave function at t = 0. So, contrary 
to the QHE where a well defined set of low energy levels are separated from the rest 
[1UI13], there is no apparent low-energy structure in the spectrum but a single "ground 
state". One could have naively expected the entanglement gap to close when reaching 
the critical point at t = 0, but this is not the case. As can be seen in Fig. [T2| the 
entanglement gap remains finite all the way from t = to t = 1 (it vanishes only at 
t = 00) We have for instance A = 1.32314 at t = 1 (exponentially fast convergence 
as a function of L) and A = 21n(7r) ^ 2.29 at t = Oj]]] A possible interpretation 
is the following: the entanglement spectrum is indeed related to the spectrum of the 

II This analytical result for A in the thermodynamic limit of the square lattice can be obtained by 
noticing that the configuration with the highest probability is | tt ■ ■ ' t) while the next configuration 
has two consecutive flipped spins Ittiit ■ ■ ■)• One can check that, for t — 0, the ratio pi/pmax of these 
two probabilities is nothing but the square of the probability for a bond located at the edge of a semi 
infinite square lattice to be occupied by a dimer. The latter probability has been computed in Ref. |20] 
and is equal to I/tt, which gives A = — ln(pi/pmax) = 21n(7r). 
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excitations that would propagate along an edge. However, in the dimer systems we 
consider, there are no gapless edge excitations, even though the bulk may be gapless for 
t = 0. 

In the thermodynamic limit, it is possible to adopt a microcanonical point of view 
where the entropy S{e) is simply related to the density of states: 

S{e) = ln(p(e)) (44) 

with 



p(e) = J]5(e-E,/L). 



(45) 



Knowing the entropy S{e) from the spectrum, the energy e(n) can be obtained as a 
function of the Renyi index n by inverting 
dS 



de 



n{e). 



(46) 



The entropy Sn is then obtained as 

5„ = ln(p(e(n))). (47) 

We conclude that, for sufficiently large L the entropy only depends of the density of 
some high energy E = L ■ e(n) in the spectrum. 






e=(E-Eo)/L 



e=(E-Eo)/L 



Figure 13. Logarithm of the density of states p associated to the entanglement 
spectrum of a half-infinite cylinder, as a function of the "energy" per site e = 
{E — Eq)/L (arbitrary units). Top: t ~ Q (square lattice), bottom: t ^ 1 (triangular 
lattice) . To display the energy range which contribute to the Von Neumann entropy 5*1 , 
the probability distribution p[e) ~ p(e) exp(— neL) is also plotted for ii = 1. System 
size: L = 28. 



The microcanonical entropy per site S{e)/L is displayed in Fig. 13 for the triangular 
and square lattice (half-infinite cylinders with L = 28). Some (finite-size) oscillations are 
visible in the triangular case, and can be interpreted as the successive energy "bands" 
corresponding to 0,2,4, ■■ ■ spin fiips in the boundary state. These oscillations will be 
smeared out in larger systems however. 
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5. Long strip geometry 

The triangular lattice can also be constructed with open boundary conditions in the 
X direction. The geometry is no longer that of a cylinder but a long strip. In such a 
situation the leading term in the entropy is still proportional to the width of the strip 
Lx = L, but the sharp corners also contribute to the sub-leading constant and it is 
not possible to extract the topological entropy for t > 0. The critical case is more 
interesting, because the first subleading correction is now a logarithm of the width. The 
later was originally predicted to be — ln(L)/4 by Fradkin and More [25] (an application 
of the Cardy-Peschel formula [H] which describes the universal logarithmic contribution 
of sharp corners to the free energy in a CFT). These terms have recently been observed 
numerically in the closely related Shannon entropy of open critical spin chains 



In Fig. 14 we show the coefficient of the ln(L) term as a function of the Renyi index n for 




Figure 14. Coefficient of tfie logaritfimic term in tfie Renyi entropy for tlie strip 
geometry, as a function of tlie Renyi parameter n. Tfiis term is extracted from a fit 
Sn = aL + b\nL + c + d/L on tlie systems sizes L — 6, L — 4, L — 2, L. Three values 
L — 14, 26, 38 are shown. The data is consistent with the CFT results. For n < ric — 1, 
the logarithmic contribution is approximately ^ —0.25 (see [IB])- For n > ric it is close 
to zero as discussed in Ref. [5^ . 

the square lattice dimer wave function with open boundary conditions. The prediction 
of Fradkin and More, — |, is verified up to n ~ 1. For larger values of n the logarithmic 
term vanishes. This is a manifestation of the boundary phase transition discussed in 
Ref. [32] • Indeed, above the compactness of the height field can no longer be ignored 
since a vertex operator cos{dh/r) (with d an integer) becomes relevant at the boundary. 
The value of d can be obtained by looking at the microscopic configuration |vax) with 
maximal probability. Contrary to the case of the XXZ chain, this configuration is non- 
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degenerate: d = 1 m the notation of Ref. [32] ■ Since the Luttinger parameter R is 
equal to 1 for the dimer problem (free fermions), the analysis of Ref. [32] immediately 
gives ric = d^/R = 1, in agreement with the present numerics. Above Uc the universal 
contribution to the entropy is that of a single "fiat" height configuration. As in the 
XXZ chain, this flat configuration does not correspond to a simple Dirichlet boundary 
condition around A in the continuum limit. Indeed, the (coarse grained) height is 
shifted by an amount 5 = ^vrr with respect to the vertical boundaries of the lattice (see 



Fig. 15). As in the XXZ chain situation, this height shift produces a logarithmic term 



which exactly compensate the logarithmic terms coming from the Cardy-Peschel angles, 
hence the absence of logarithm in the Renyi entropy when n > ric = 1. 



If 



max' 



B 



-1 




Figure 15. Configuration |imax) with the maximal probabihty on the square lattice 
and a compatible dimer covering of the rectangular region A. The microscopic heights 
are indicated in units of ^Trr. When turning clockwise around a site of the even (resp. 
odd) sublattice the height changes by +1 (resp. -1) when crossing an empty bond, and 
changes by -3 (resp +3) when crossing a dimer. The lower horizontal boundary of A has 
a coarse-grained height which is "flat", with an average height equal to |(0 + 1) — \ 
(red). The vertical boundaries have a coarse-grained height equal to ^(1 + 2) = | 
(green). In the continuum limit there is an height shift 5 = at each corner of A. 



6. Kitaev-Preskill construction 

As discussed in Sec. |4] the cylinder geometry allows to extract the subleading entropy 
term in a rather straightforward way, by a simple fit of Sn{L) on (at least) two system 
sizes. However, the original proposals [3 E] were to extract the topological entanglement 
entropy from a single and large planar system. There, the subsystems on which the 
entanglement entropy are computed cannot have a straight boundary and necessarily 
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have comers, etc. These comers (as well as the curvature) also contribute to the 
entanglement entropy by a (non-universal) amount of order one and therefore need 
to be subtracted. The subtraction scheme proposed by Kitaev and Preskill [6] is based 
on the following combination of entropies (see Fig. 16). 

r-topo _ n(ABC) _ q(AB) _ r.(BC) _ qiAC) , r-CA) , q(B) , q{C) / ^a^ 

'-'n — '-'n '-'n ^ '-'n ^ '-'n ^ \^°) 

The first numerical implementation of this subtraction ideas was done in a the RK 
dimer wave function at t = 1 and n = 1 [25]. Some other recent works investigated the 
n = 2 case using quantum Monte Carlo on a Bose-Hubbard model and variational 
quantum Monte Carlo on projected spin liquid wave-functions [16]. Here we extend 
the results of Ref. [2S] on dimer RK wave functions for several values of t, n, and with 
with finite areas A, B and C embedded in a infinite plane. The results are shown in 
Fig. 17 Provided t is not too small [i.e. the dimer-dimer correlation length is not too 
large), the Kitaev-Preskill construction gives an entropy constant equal to — ln(2) with 
high precision, as expected. Still, for the same numerical effort (boundary length), the 
convergence turns out to be slower than with the cylinder geometry (see Fig. [s]). 



Figure 16. Geometries required for the computation of sl!^^'^\ sli^^^ and S'i'^' at 
Radius p = 4.5. They have Ni, — 30, 29 and 19 boundary sites (in red) respectively. 



The Eq. 48 was originally designed to probe massive wave-function, but it is also 
natural to consider the limit t — ?■ where the wave function becomes critical (and 
restricting to n < Hc for simplicity). 

Each term in Eq. [48] corresponds to a subsystem Q = ABC, AB, ■ ■ ■ which is 
topologically equivalent to a disk, but possibly with some sharp corners. For each 
such subsystem, we wish to use a formula derived in Ref. [32] : 



n 



In 



\ UK 



(49) 



where ^ is a free-field partition function on the whole system, and is the partition 
function with Dirichlet boundary condition imposed at the boundary of Vt (thus 
disconnecting VL and 0). k is the bare stiffness and the first term should be evaluated 
with a modified stiffness k' = 

By construction, the non-universal contributions proportional to the boundary 
length will drop out of the KP combination. Next, we consider the logarithmically 



% This formula was originally derived in the case in the case where (7 is a half infinite cylinder, but the 
argument in fact applies to the present geometries as well. 
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Figure 17. Top left : critical case t ^ Q. Top right : t = 0.3. Bottom left t = 0.7. 
Bottom right t = 1. In each case, — iS'*°P°/ln2 is shown, for n = 0.5,0.8, 1, 1.5,2 as a 
function of the radius p. 



divergent terms which come from the sharp corner contributions to the free energies. 
Each corner with interior angle a gives a contribution F{a) = ^^(f - ^) ln(L//o) to the 
free energy, where L is the typical scale of the boundary and /q some microscopic cut 
off IS]. To apply the Eq. 49 what needs to be computed is the free energy difference 
between that of the whole system, and that where VL and VL have been disconnected 
(Dirichlet boundary condition). So, in the disconnected term, a sharp corner of angle a 
in Vt will also contribute as a sharp corner of angle 27r — a (in Vt). The contribution to Sn 
is thus 5Sn = F{a) + F{2tt — a) = ^(2 — ^ — 2^^^) ln((L//o), which is by construction 
symmetric under the exchange a O 27r — a. Then it is easy to check that in the spatial 
decomposition implied by 48, each angle appearing in some +Sn will cancel out with 
another one (with the same angle or its complement) in —Sq/. 

However, as already mentioned in Ref. |15], this is only true for the leading 
(logarithmically divergent) part, because there is no simple reason why the microscopic 
length scales Iq should all be the same. We thus expect some constant (non-divergent) 
and non- universal contribution to the entropy when t = 0. 

Refs. mentioned that the entanglement entropy of a disk Q of radius R 

embedded in a larger disk Q of radius L could have a (very slowly) diverging term 
~ ln(ln(L/i?)) for a critical RK wave function. However, in the lattice (dimer) version 
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of the RK state we consider, it is easy to show that the entropy must be finite when 
L — oo while keeping R fixed. The argument is as follows: the (Von Neumann) 
entropy 5*1 of a subsystem can be expressed using the probabilities pi of its boundary 
configurations: 

N 

Sl = -Y,V^HVi) (50) 

i=l 

where N is the number of possible microscopic configurations at the boundary of Vt. 
If the boundary has a finite length R, N must be finite with InA^ ~ i?. As a 
consequence, since the entropy is bounded by InA^, we have Si < R. In other words, 
the entanglement entropy cannot exceed the boundary law for RK states. This bound 
does not involve the size of the outer system Q, and none of the entropies appearing in 
Eq. 48 can diverge when taking the outer system to its thermodynamic limit. Why the 
argument of Ref. [2S] does not apply to this quantity in lattice RK states is however 
unclear to us. But in any case S'topo(-R) cannot diverge when taking L — )• oo at fixed R, 
whatever the lattice RK state provided it has a finite number of states per site. This is 
indeed confirmed by our numerical estimations of 5'topo(-R) which are performed directly 
in the thermodynamic limit L = oo and which gives finite values for finite values R. 
Although the system sizes {R) are too small to observe the true large-i? behavior for 
t = (square lattice), the argument above concerning the corner contribution indicate 
that it is very likely a non-universal number. 



7. Summary and conclusions 

Thanks to some extensive use of the Pfaffian solution of the classical (2d) dimer model, 
we have performed exact calculations of the entanglement entropy and entanglement 
spectra of some dimer RK states on large subsystems. Using the cylinder and the 
Kitaev-Preskill geometries we recovered the topological entanglement entropy of the 
Z2 phase, — ln(2), with high accuracy. As expected, this value not only holds for the 
triangular lattice RK wave-function, but is in fact independent of the fugacity t > 0. 
We also analyzed the scaling close to the critical point at t = 0, as well as the behavior 
for large values of the Renyi index n. In particular, we proved for n — t- 00 that the 
sub-leading entropy constant is — ln(2). Thanks to its translation-invariant boundary, 
the cylinder geometry gives smaller finite-size effects and therefore a much more precise 
estimation of the topological entanglement entropy than the KP setup (for a given length 
of the subsystem boundary). For this reason, it may be preferred in future numerical 
studies (exact diagonalization or quantum Monte Carlo) looking for topological ground 
states in realistic lattice models. 

The entanglement spectra were also computed in the cylinder geometry, and the 
presence of a unique ground-state and a finite gap (whatever the fugacity) showed that 
for these states, contrary to naive expectations, the topological (or critical) nature of the 
phase is not apparent in the low-energy part of the entanglement spectrum. Simpler Z2 
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wave functions such as that of the Toric Code [13] (or that of Ref. [H]) do not allow to 
learn much about the structure of the entanglement spectrum. Indeed, in those states 
with vanishing correlation length all the non-zero eigenvalues of the reduced density 
matrix are exactly degenerate (no n dependence of the Renyi entropy). From this 
point of view, the dimer states we consider offer an interesting compromise between 
the possibility to do exact calculations on large systems and a non-trivial entanglement 
spectrum. Extending these calculations to other states with richer topological structure, 
like string-nets wave functions ^47J, could be a promising direction of research. 
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Appendix A. Green function elements for an infinite cylinder 

Appendix A.l. Diagonalization of the Kasteleyn matrix 

We wish to diagonalize the Kasteleyn matrix by Fourier transform for — )■ oo. To do 
so we must distinguish between two sublattices (see Fig. Isj): 



£o = {(2xux + yuy) I < a; < LJ2, < y < L J (A.l) 

£i = {(2x + l)ux + yuy I < X < L^./2, Q<y <Ly} (A.2) 

We denote by = L^Ly the number of sites. Then we define a new basis 

IM)^^E^-'W (A.3) 

V / roeCo 

|k,l) = — ^ V e-^'^-'-Mri) (A.4) 

V / ri6£i 

The Kasteleyn matrix satisfies antiperiodic boundary conditions in the x— direction, 
and since Ly — )• oo, we can also assume antiperiodic boundary conditions in the y— 
direction. The appropriate wave-vectors are the k = fc^^Ux + kyUy with 

(2j + l)7r 



kx G Kx 



ky G Ky 



Lx 
(2j + l)7r 



Ly 



j = 0,...,L,/2-l| (A.5) 
j = 0,...,L,-ll (A.6) 



In the new basis, the Kasteleyn matrix takes the following simple form 

^ ilA - ( 2ismky 2i sin kx + 2tcos{kx + ky) \ 

'^^ I 2i sin kx — 2t cos{kx + ky) —2i sin ky j ' 

and can easily be inverted 

^_i/j^N 1 ( —2ismky —2i sin kx — 2t cos{kx + ky) | g\ 

~ det [JCaM] \ -^isinkx + 2tcos{kx + ky) 2isinky J 

with 

det [Ka/sO^)] = A sin^ kx + 4: sin^ ky + At"^ cos'^{kx + ky). (A. 9) 

For two sites r = xUx + yUy and r' = x'ux + y'uy respectively in sublattices a and (3, 
the Green function element is 

1C;1, = ^ V e-'^^("'-") / dky IC'l (k)e-''=-(^'-^) (A. 10) 

T^Lx Jo 

In this equation, the integral on dky can in principle be done explicitly for any y' — y, 
as will be shown in the next subsection. To compute the entanglement entropy in the 
cylinder geometry \y' — y\ doesn't however need to be greater than 2, whereas it can 
attain 3 in the strip geometry. 
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Appendix A. 2. Green function elements 

The computation of Green functions element requires the evaluation of integrals of the 
form 

C{k)= r COS{pky) 

Jo Asin^kx + Asin^ky + At^cos-^ikx + ky) ^ ^ ' ' 

^pl/c.j - 4smHx + Asm'ky + At^cos\k, + ky) ^ ' 

with p an even integer (otherwise the integrals are simply zero by symmetry). Both 
integrands are tt— periodic and following Bioche's rules we can make the change in 
variables u = tan ky. We get 

C(k)=' r T,[(l + .rV^] du 

2J_^ m2[1 + (1 +t2)sin2A;J -Mt2sin(2fc^) + sin2A;^ + t2cos2A;^ ^ ' ^ 

c., r u{l + u-)-y-Up.,[{l + u-)-^l-]du 

2y_^ u2[l + (l+t2)sin2fc,] -Mt2sin(2A;,) + sin2fc, + t2cos2fc, ^' ^ 

where Tp{x) and are the Chebyshev polynomials of the first and second kind 

respectively: 

Tp{cos9) = cospe (A. 15) 

Up.^{cos9) = ^-^ (A. 16) 

smfc* 

For p even Tp{—x) = Tp{x) and f/p_i(— x) = —Up^i{x). Therefore, both integrands in 
Eq. A. 13 and A. 14 are rational functions of u, as should be. Cp and Sp can then be 
calculated by residue. Closing the contour by a big circle in the upper-half plane, two 
poles will contribute to the integral. The first pole is at 

^2 sin k^ cos k^ + i ^/f^ + sin^ k^ + sin^ k^ , . _s 

u = ; (A. 17) 

1 + (1 + t2) sin^ kx ^ ' 

and is of order 1. The second one at -u = i is there if p 7^ and is of order p/2. 

Although the residue calculation for any even p is in principle straightforward, the 

procedure becomes more and more cumbersome when p gets bigger. Only for p = do 

we get a simple (known[16]) result: 

Co{k.x) = . = (A.18) 

V t"' + sm k^ + sm k^ 

From these we can get access to all the Green functions elements. The simplest are 

along the same horizontal line, and only require the knowledge of Cq : 

^2lu. =0 (A. 19) 

r-i _ 1 ^ sinfc^sin(2£+ l)fc^ 

^{2£+i)u, - r / ,2 I -2, — I ■ 4, ^^-^^^ 

J^x \/t^ + sm kx + sm kx 

For the cylinder geometry, the knowledge of Cq, C2 and 5*2 is sufficient. For the strip 
geometry, also C4 and 5*4 are needed. To compute the entanglement entropy in the 
Kitaev-Preskill geometry, it is easier to evaluate the double integral (L^. — )■ 00) in 
Eq. A. 10 numerically. 
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Appendix B. Closed-form formula for 5'„-oo in the cylinder geometry 

As explained in the text, the maximum probability corresponds to a simple configuration 



with all boundary spins up. Then, a natural way to proceed would be to use Eg. [27] and 
try to evaluate the resulting determinant. This method is most certainly viable, but we 
will follow another path. In the dimer language, the probability we are looking for is 
given by 



lim 



[Zcy\{Lx, Ly/2)] 



7 (r r ) ^ (^-^^ 

where Zcyi{Lx, h) counts the number of dimer coverings on a finite cylinder of 
circumference and height h. Despite the loss of translational invariance in the y— 
direction, Zcyi can still be evaluated in closed form, as is shown in Appendix B.l From 
this Pmax can easily be calculated, see Appendix B.2 



Appendix B.l. Dimer coverings on a finite cylinder 

Let Zcyi be the partition we are looking for. Using (skew) translational invariance along 
the X-axis, one gets (recall = {{2m — 1)t{ /L^ , 1 < ?ti < Lj./2}): 



Zcyi(L^., Lyf = Y\ det 



^<iJ<'^Ly 



(B.2) 



In other word, the Kasteleyn matrix is block-diagonal with Lx/2 blocks of size 2Ly. 
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Although it is not easy to diagonahze K,^^^ , its determinant can be exactly evaluated 
using the perturbation trick. To do so, we introduce 
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-tx 1 

-1 -tx 
\ t, 1 

This amounts to putting antiperiodic boundary condition along the y— axis for the total 
Kasteleyn matrix. /Cq^^ is block skew circulant, and it can be diagonalized in Fourier 
space. In particular its determinant can be easily evaluated : 

det4") = H A{k„ky) 



(B.5) 



A{k^,ky) = A sin^ k^ + A sin^ ky + At^cos^{k^ + ky), (B.6) 



where Ky — {{2m — l)n/Ly ,1 < m < Ly}. This allows to express det/C*^^) as 



det/C(^) 



det/C 



{x) 



det 1 + 



(x) 



detM 



(x) 



(B.7) 



)(^(.^) — is a matrix with only 8 non-zero elements, and using elementary row-column 
manipulations, the determinant can be reduced to a 4 x 4: 

/ z —a w —ib \ 
—a z ib —w 

—w ib z a 
—ib w a z 

After some algebra, we get the following formulae for the coefficients : 



M, 



(x) 



\ 



{z, w, a,b) eCxCxRx 



I 



1 _^ 2 ^ sin^ k^ i [sin(2A;j,) - sin(2^^ 2ky)\ 



2t 



It X - 
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COS k 



A.[kx, ky) 



A.[kx, ky) 



2it sin kx e ^'^^ 



y h 
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A.[kx, ky) 

sin kx 
A{kx, ky) 



(B.9) 
(B.IO) 
(B.ll) 
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The number of dimer coverings on the triangular lattice with cylindrical boundary 
conditions is then given by: 

Zeyi(L,,L,) = n ] det (Mf ^) X \{^{k,,ky) \ (B.13) 

Evaluating the determinant, we finally get the following closed formula for the partition 
function 



(B.14) 



with 



{t^ + sin^ + sin^ k^ dik^)"^ + (sin^ k^ — t cos kj.)d{kx) + 1/4 + 
2 ^ 1 

2 ^ sin(2A;j^) - sin(2A;^. + 2ky) 



A(fca;, ky) 



(B.15) 



Appendix B.2. Exact formula for Sn=oo 

The maximum probability is in the thermodynamic limit given by 

[Zcyl{Ly/2, Lx)f' 

Pmax= lim 



Ly^OO Zcyl{Ly, La 



(B.16) 
(B.17) 



Eq. B.17 follows from Eq. B.16 using Euler-Maclaurin's formula on the ratio of terms 
involving A{kx, ky), coming from Eq. B.13, Using Eq. A. 18, we also have 



1 



lim d{kx) = — -. 

^y^"^ 2Jt^ + sin^ k^ + sin^ k^ 



(B.18) 



while lim e{kx) = because the integrand has a symmetry center solution of sin(2fcj^) 

Ly—^OO 

sin(2A;a- + 2/cjy). In the end we obtain 



l<m<L/2 



5*00 = -Inpn 



E, / 1 1 sin^ kx — t cos k^ 
In - + 



_ (2m-l)7r 



(B.19) 



Appendix B.3. Asymptotic expansion 

At t = 0, the subleading constant in the L — )■ oo asymptotic expansion just follows from 
the Euler-Maclaurin formula. We find 



Soo(t = 0) =0. 



(B.20) 
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Some additional care must be taken in the case t > 0. The function 

, /I 1 sin^ A; — t cos A; \ ,^ 
f{k) = - In - + -—====== (B.21) 

\2 2 ^lYi^ k + siiC k J 

actually diverges as f{k) ~ — 21n/c - independent on t - when — )■ 0. The asymptotics 

can be obtained by applying the Euler-Maclaurin formula on [f{k) + 21nfc] while 

applying Stirling's formula on the remaining "linearized" term — 2 In k. Doing so 

we finally obtain the topological term 

> 0) = -ln2. (B.22) 

Only the linearized term actually contributes to the constant. Indeed, it is universal 
and shouldn't be affected by the short-distance {i.e high momentum k) details of the 
model. 
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